Evolutionary phase space in driven elliptical billiards 
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We perform the first long-time exploration of the classical dynamics of a driven billiard with a four 
dimensional phase space. With increasing velocity of the ensemble we observe an evolution from a 
large chaotic sea with stickiness due to regular islands to thin chaotic channels with diffusive motion 
leading to Fermi acceleration. As a surprising consequence, we encounter a crossover, which is not 
parameter induced but rather occurs dynamically, from amplitude dependent tunable subdiffusion 
to universal normal diffusion in momentum space. In the high velocity case we observe particle 
focusing in phase space. 

PACS numbers: 05.45.Ac,05.45.Pq 



Billiards are a widespread paradigm to study complex 
dynamics in many areas of physics. Two dimensional 
systems are especially appealing, since unlike one dimen- 
sional ones they possess already a rich phase space struc- 
ture, which can range from integrable over mixed to fully 
chaotic Billiard systems can be realized experimen- 
tally in different ways, such as using semi-conductor het- 
erostuctures [H, quantum dots atom optical [H and 
(superconducting) microwave setups 0, 0] . Very recent 
applications of their dynamics are the investigation of 
directed emission from laser-microcavities @, H, Q and 
the design of improved thermoelectric efficiency, where 
billiards are used to tailor the desired microscopic prop- 
erties [io| . 

A natural generalization are time-dependent billiards, 
where the boundary is driven, for example harmonically. 
This allows the study of non-equilibrium processes such 
as Fermi acceleration (FA) [ll|, [l2|, [l3[ . There are only 
few studies of such two dimensional (2D) driven billiards 
HI, 15, 2, 17, 3, HI, [2(j, since the long-time numerical 
iteration of the corresponding four dimensional (4D) im- 
plicit mapping is computationally challenging and in par- 
ticular the resulting phase space is hard to analyze due to 
its high dimensionality and unboundedness. So far this 
prevented any investigation addressing the long-time be- 
havior of corresponding dynamical quantities. Key ques- 
tions are whether there will be FA pjj, i.e., unbounded 
energy gain of an ensemble exposed to time-dependent 
external potentials and more specifically what laws this 
FA will follow depending on the geometry and the driv- 
ing of the billiard. In one dimension, the existence of FA 
depends exclusively on the driving law [2l[, in 2D the 
situation is much more complicated, since the dynam- 
ics of the corresponding static systems can be integrable, 
mixed or chaotic. Recently the existence of FA in the 
harmonically driven elliptical billiard was demonstrated 
(20| . augmenting the "LRA-conjecture" which says that 
only time-dependent systems with a mixed or chaotic 



static counterpart exhibit FA [15]. However, none of the 
mentioned works analyzed the full 4D phase space and 
investigated how its structure determines the diffusion 
in momentum space, especially with respect to the long 
time behavior of the diffusion and acceleration process. 

In this Letter, we present, for the first time, a detailed 
exploration of the long-term dynamics of a driven bil- 
liard with a 4D phase space. It is shown that the phase 
space geometry evolves, with increasing velocity of the 
propagating ensemble, from a large chaotic sea with rich 
stickiness structures due to regular islands, to large dom- 
inating regular domains with thin chaotic acceleration 
channels. This change in the phase space structure leads 
to a crossover from amplitude dependent tunable subdif- 
fusion to universal normal diffusion in momentum space. 
Due to its acceleration with increasing number of colli- 
sions the particle ensemble propagates effectively in an 
evolutionary phase space which we study in depth. 

The geometry of the time-dependent elliptical billiard 
is given by 

/ x(t) \ _ / a(t) cos _ ( (Aq + C sm(ujt)) cos <j) 
\y\t) ) ~ V Kt)sm<j> /V l B o + Csin(ut))sin0 

(1) 

where (x(t),y(t)) is a point on the boundary, <\> is the 
corresponding elliptic polar angle, C > is the driv- 
ing amplitude, uo is the frequency of oscillation and Aq, 
Bo are the equilibrium values of the long and the short 
half-diameter, respectively. We fix uj = 1, Ao = 2 and 
Bo = 1 without loss of generality (arbitrary units). Due 
to ballistic motion in between collisions, the dynamics 
can be described by an implicit AD map using the 
variables (£ n , n , a n , v n ), where v n = \v n \ is the modulus 
of the particles velocity, £ n = cot mod 2tt is the phase of 
the boundary oscillation and a n is the angle between v n 
and the tangent on the boundary at the n-th collision. 
Note that the corresponding static ellipse is integrable 
[l|, since besides the energy the product of the angular 
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where e is the eccentricity, is conserved. The static phase 
space is globally divided by the separatrix {F = 0), con- 
necting the two hyperbolic fixed points, into rotators 
(F > 0) and librators (F < 0). Two elliptic fixed 



points are located at the minima F n 



7(i -£ 2 ). 



Periodic orbits with short periods are located around the 
separatrix, whereas the existence of period orbits in the 
librator region close to the elliptic fixed points requires 
comparably long periods [14]. This will later be impor- 
tant for the understanding of the distribution of periodic 
orbits in the driven case. 

In ref. [20[ we showed that the driven ellipse exhibits 
Fermi acceleration (FA). Here we present a comprehen- 
sive analysis of the AD phase space, specifically we ex- 
plore how its composition changes with increasing v and 
demonstrate how this affects the long-time evolution of 
FA in the system. We first analyze the periodic orbits and 
their stability and consider subsequently the phase space 
density p of a typical ensemble in detail, by exploring its 
dependence on the number of collisions n as well as the 
velocity v. Due to the occurrence of FA, the phase space 
is open with respect to v. Unlike in the static system, a 
is not restricted to [0, 7r], but can now range from to 
2tt [3]. Nevertheless, the total phase space is not just 
given by the 4D cuboid [0, 2tt] x [0, 2tt] x [0, 2tt] x [0, oo] 
spanned by £, (/>, a and v, but possesses a more compli- 
cated topology, since e.g. for £ G [tt/2, 3tt/2] (the ellipse 
is contracting) a is restricted to a G [0, tt] [221 ] . 

Finding (unstable) periodic orbits in such a high di- 
mensional phase space is an intricate task and becomes 
more difficult with increasing velocity of the orbits, since 
the typical period grows linearly with v. We use a vari- 
ant of the method developed in Since the mapping 
of the driven ellipse is 4D and area-preserving, there are 
three possible types of eigenvalues of the corresponding 
real monodromy matrix: (i) Four complex eigenvalues 
(elliptic), (ii) two real and two complex eigenvalues (lox- 
odromic) and (iii) four real eigenvalues (hyperbolic). The 
first case corresponds to stable [25 L w hereas the latter 
two represent unstable fixed points [2l(. Performing ex- 
tensive numerical computations we detected the periodic 
orbits up to period 100. The density p(F) of the peri- 
odic orbits as a function of a) is shown in Fig. [H 
By using F, we can effectively map the (j) x a plane onto 
a single dimension (see the insets of Fig. [T] where the 
F = const, isolines are shown in black). Fig. 1 shows 
periodic orbits with v\ < v < V2 for (a) v\ = 0, V2 = 1 
and (b) v\ = 30, V2 = 40. For small values of v (see Fig. 
□JO, the periodic orbits up to period 100 cover the whole 
available F-range, and thus the whole <p x a plane, ex- 
cept for a narrow region F « around the separatrix of 
the static system. This has to be expected, since the col- 
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FIG. 1: (Color online) Density of periodic orbits up to period 
100 as a function of F with v\ < v < V2 and v± = 0, v<i = 1 
in a) and v\ — 30, V2 = 40 in b). The F(<f>,a) ~ region 
is depleted for low velocities (a), and is populated predomi- 
nantly for higher velocities (b). Exemplary, the distribution 
of periodic orbits is shown in the insets in the <j) x a plane, 
where the black curves are the F — const, isolines. 



lisional dynamics with the oscillating walls deviates the 
most for low v from the dynamics of the corresponding 
static system. The region of the separatrix is associated 
with the highest instability and in particular the stable 
periodic orbits existing in the static case around F w 
have been destroyed. The opposite scenario can be seen 
for high values of v, see Fig. [TJd. Here we expect the 
fingerprints of the emergence of a quasistatic regime for 
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FIG. 2: (Color online) Collision resolved phase space den- 
sity pn 1: n 2 (F,^) (log colormap), m — 10, 712 = 10 2 (a) and 
n\ = 4 • 10 6 , 712 = 10 7 (b). p ni) n 2 (i ? , £) is largest around 
F w —2.5 for small n (a), whereas the F{<j>,a) ~ region is 
predominantly populated for high collision numbers (b). In 
the insets, the corresponding p ni ,n 2 (0? a ) are shown. 
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the billiard: The particles are fast enough to experience 
many collisions with the oscillating walls within a narrow 
interval of the phase of oscillation. The relative momen- 
tum change obeys Sv/v <C 1 at each collision and the 
angle of incidence is almost the angle of reflection. As 
a result they trace orbits of the static system. As can 
be seen in Fig. [TJd, the region F « possesses the high- 
est density of period orbits, in particular unstable ones 
(hyperbolic and loxodromic), whereas the librator region 
around the two elliptic fixed points of the static system 
is depleted. 

The consequences of the above described changes in the 
distributions of the periodic orbits with increasing v be- 
come immediately clear when studying the evolution (in 
terms of the number of collisions) of the phase space den- 
sity p ni ,n 2 (F, = U2 _^ 1 + 1 Enlnx Iv P( F > & V ' U ) dv ° f an 

ensemble 

1 n p 

p(F,t,v,n) = —^5(Z-g)5(F-F(W,a?))5(v-v?^ 
p i=i 

(3) 

where the indices stand for the zth particle at collision 
number n, respectively. An ensemble of N p — 10 3 par- 
ticles, with vo = 0.1 and ^0,^0,^0 distributed uniformly 
and randomly (£0, 4>o £ [0, 27r], ao G [0, 7r]), is iterated 
for 10 T collisions. The visiting probability is shown in a 
log scale colormap in Fig. [2j In the beginning (n\ = 10, 
n 2 — 100) of the evolution, the region around the ellip- 
tic fixed points of the static system possesses the highest 
visiting probability, see Fig. [2^ and especially the inset, 
where the corresponding p ni ,n 2 (0? °0 is shown. For large 
collisions numbers (n\ = 4 • 10 6 , n 2 = 10 T ), see Fig. (2}d 
and the inset, the ensemble is located predominantly in 
a region around the separatrix (F((j), a) = 0) of the cor- 
responding static system. In the insets, the half-space 
a > 7r is not shown, since p(0, a) ~ there. 
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FIG. 3: (Color online) Velocity resolved phase space density 
Pv!,v 2 (F,£) (log colormap). For low velocities (a), the ensem- 
ble covers the whole Fx^ plane, for higher velocities (b), it is 
exclusively located around the F(<p,a) ~ region. The insets 
show the corresponding p VljV2 ((f), a) . 



The effect of the population of the region around 
the separatrix becomes even more obvious when con- 
sidering the phase space visiting probability velocity 
resolved rather than collision resolved, p VlyV2 (F^) = 
n7 Ivi En P(F, £, v, n)dv. All collisions between m = 1 
and n 2 = 10 7 with v\ < v < v 2 are projected onto the 
F x £ space. For low velocities, v\ = 0, v 2 = 1 (see 
Fig. [3fci), the complete F x £ space and thus the com- 
plete (j) x a plane (see the inset where the corresponding 
p Vl iV2 ((p,a) is shown) is populated. For higher velocities, 
v\ = 30, v 2 = 40 (see Fig. [3}d and the inset), a narrow, 
sharply bounded region around F(<p : a) = is now ex- 
clusively populated, i.e., outside this region no collisional 
events occur. This kind of population of phase space is 
generic and does not depend on the initial ensemble, as 
long as the initial velocity vq of the ensemble is suffi- 
ciently small. Unlike in dissipative systems the above 
squeezing of the ensemble is not caused by the existence 
of attractors but is due to a mechanism which will be 
explored in the following. 

For low v we encounter a mixed phase space, where 
an infinite hierarchy of regular islands is embedded into 
a large chaotic sea. With increasing v, large regular re- 
gions emanate (we confirmed the existence of these re- 
gions among others by calculating the corresponding Lya- 
punov exponents and determined their extension) from 
the position of the elliptic fixed points of the static sys- 
tem. Their extension in the <p x a subspace grows rapidly 
with increasing v until v « 10. From there on the reg- 
ular regions capture a substantial amount of the avail- 
able phase space volume, except a small region around 
F((j), a) w 0. If we now increase v further, the growth of 
these regular regions is reduced significantly, leaving the 
channel F(<fr,a) « open for diffusive processes. The 
precondition to use sufficiently low vo in order to observe 
FA becomes obvious now: It ensures that particles start 
inside the chaotic sea which is connected to the channels 
F(4>, a) w that eventually allow for the acceleration 
process. 

The evolution of the phase space density can then be 
understood as follows. The pronounced localization of 
the ensemble with respect to the F x £ and correspond- 
ingly to the (j) x a space in the v resolved picture (Fig. [3|) 
is smeared out when considering the collision resolved 
density (Fig. [2]). This happens because the increase of 
the ensemble average (v) with increasing n comes along 
with a certain distribution p(v) (a Maxwell-Boltzmann 
like distribution [22[), i.e., the n resolved distribution is 
a weighted superposition of the v resolved Pv lt v 2 (FiQi 
where the weights are given by p(v). 

The evolution of the composition of phase space with 
increasing v has amazing consequences on the mean ve- 
locity (\v\) of the ensemble, in particular on its long time 
behavior. The rich phase space structure, associated with 
many spacious regular islands embedded into the low v 
chaotic sea, causes stickiness of the trajectories. As a 
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FIG. 4: (Color online) Ensemble averaged velocity (\v\)(n) ~ 
n« c > as a function of n (collisions) for two different ampli- 
tudes C . Note the different slopes (3 between ni = 3 • 10 4 — 10 8 
and n 2 = 10 8 -10 n : in the interval m /?(0.1) = 0.28, /?(0.5) = 
0.40 (subdiffusion), in the interval n 2 (3(0.1) = (3(0.5) = 0.5 
(normal diffusion). 



result the dynamics is intermittent, exhibiting both lam- 
inar phases of oscillations as well as chaotic phases. The 
chaotic phases of the trajectories lead to stochastic fluc- 
tuations of v and thus to a net increase of the mean 
velocity (\v\) of the ensemble. This intermittent inter- 
play of laminar and turbulent phases can be effecti vely 
described through continuous time random walks [24[ 
predicting the appearance of subdiffusion, which is ex- 
actly observed here. With increasing (\v\), the ensemble 
is more and more squeezed into the acceleration chan- 
nels around F = 0. A detailed analysis based on the 
periodic orbits shows that there are much fewer and at 
the same time smaller regular islands in the acceleration 
channels for higher compared to low velocities, see also 
Fig. [U This leads to a significantly reduced stickiness, 
i.e., fewer and shorter laminar phases once the ensem- 
ble is predominantly located in the narrow acceleration 
channels. Consequently, the diffusion in v space becomes 
then normal. The crossover at around n c ~ 10 8 from 
subdiffusion to normal diffusion can be seen in Fig. [H 
for two different values of the amplitude (C = 0.1 and 
C = 0.5). Such a crossover occurs in other driving modes 
of the ellipse as well, e.g. in the quadrupole mode, where 



a(t) = ao + Csinujt and b(t) = bo — Csinutt [22[. In 
Fig. [4j the average velocity grows according to a power 
law (\v\)(n) ~ between n « 5 • 10 4 and n « 10 8 . 

In this region, the diffusion exponent /3(C) is amplitude 
dependent (/?(0.1) = 0.28, /?(0.5) = 0.40) and always 
smaller than 0.5 (subdiffusion). The dependence of /3(C) 
in this regime can also be traced back to the phase space 
properties, e.g. there are more islands causing stickiness 
for smaller than for larger amplitudes [22]. For n > 10 8 
(\v\)(n) ~ n 1 / 2 , i.e., the exponent obeys (3 = 1/2 inde- 
pendent of the amplitude. Note however that an ensem- 
ble of particles has to be propagated for n ^> n c ~ 10 8 



(until n = 10 11 in our case), which is extremely demand- 
ing and requires large scale computations, to detect this 
crossover. We checked the convergence of the ensemble 
with 160 particles shown in Fig. [H with a reference en- 
semble consisting of 10 4 particles and shorter iteration 
times and found perfect agreement. 

Our analysis of the 4D phase space of the driven el- 
liptical billiard shows how its composition dynamically 
changes with increasing velocity v. For high v, thin 
acceleration channels of diffusive motion evolve from a 
spacious chaotic sea located at low v. This induces a 
crossover from sub to normal diffusion in momentum 
space. The acceleration channels are present in particu- 
lar in the high v quasistatic limit. We expect this to be a 
generic feature of many driven dynamical systems which 
show Fermi acceleration. With increasing v, the ensem- 
ble will be focused on phase space regions where the cor- 
responding static system possesses the highest degree of 
instability, or even chaoticity in the case of non-integrable 
billiards. In this sense driven billiards, as shown in the 
case of elliptical geometry, may operate not only as an 
accelerator but also as a phase space collimator. 
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